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We have carried out extensive Monte Carlo simulations of the fusion of tense apposed bilayers 
formed by amphiphilic molecules within the framework of a coarse grained lattice model. The fusion 
pathway differs from the usual stalk mechanism. Stalks do form between the apposed bilayers, but 
rather than expand radially to form an axial-symmetric hemifusion diaphragm of the trans leaves 
of both bilayers, they promote in their vicinity the nucleation of small holes in the bilayers. Two 
subsequent paths are observed: (i) The stalk encircles a hole in one bilayer creating a diaphragm 
comprised of both leaves of the other intact bilayer, and which ruptures to complete the fusion pore, 
(ii) Before the stalk can encircle a hole in one bilayer, a second hole forms in the other bilayer, and 
the stalk aligns and encircles them both to complete the fusion pore. Both pathways give rise to 
mixing between the cis and trans leaves of the bilayer and allow for transient leakage. 



I. INTRODUCTION 

Although membrane fusion is a fundamental biological 
process of importance in fertilization, synaptic release, 
intracellular traffic, and viral infection, its basic mecha- 
nism is not well understood. Much of the literature has 
focused on fusion proteins whose function is, inter alia, 
to overcome the energetic cost of bringing the bilayers to 
be fused to within a small distance of one another, a step 
which places the membranes under tension 0|. There 
is accumulating evidence, however, that the subsequent 
stages in the fusion pathway, the interruption of the in- 
tegrity of the bilayers and the molecular rearrangements 
that lead to the formation of the fusion pore itself, are 
essentially lipidic in nature . A consequence of 

this view is that the fusion process can be studied, both 
experimentally and theoretically, utilizing simple model 
membrane systems. Knowledge of the fusion mechanism 
in these simpler systems would illuminate additional roles 
that the proteins need to play in biological fusion. 

The theoretical treatment of membrane fusion has, al- 
most without exception, been restricted to phenomeno- 
logical models which decribe the bilayer, not in terms 
of the microscopic architecture of its components, but 
rather in terms of the macroscopic elastic properties of 
its monolayers. The common assumption is that these 
elastic moduli are uniform and independent of membrane 
deformations Although attractive mathematically, 

this approach has its limitations. For instance, it is not 
clear whether the expansion of the membrane free energy 
to second order in deformations is sufficient to describe 
the highly curved intermediate structures, which may be 
involved in fusion. Additional approximations must be 
introduced to calculate the properties of junctions of bi- 
layers, which are not well described by simple bending 
deformations. The energy of these structures has proven 
to be particularly sensitive to the approximation used in 



their description |0, |23, . Importantly, applica- 
tion of these approaches requires one to assume a par- 
ticular fusion pathway. The only pathways considered 
to date has been limited to variants of one hypothesis 
0, 0, nil • One starts with two bilayers in close apposi- 
tion. Lipids in the facing, or cis, layers rearrange locally 
and bridge the aqueous gap between the bilayers. This 
results in the formation of an axially symmetric stalk. In 
most versions, the stalk then expands radially and the cis 
layers recede. The trans layers make contact and produce 
an axially symmetric hemifusion diaphragm. Nucleation 
of a hole in this diaphragm completes the formation of an 
axially symmetric fusion pore. Because of the evolution 
of the stalk into a hemifusion diaphragm in this model, 
we shall refer to it as the "hemifusion mechanism." Be- 
cause only variants of the hemifusion mechanism have 
been examined, and because the theory is phenomeno- 
logical, one does not know a priori in what systems this 
pathway may be the most favored, or under what con- 
ditions. Some insight is gained by comparison with ex- 
periment which shows this hypothesized mechanism to 
be consistent with a wide range of experimental observa- 
tions of biological lipids 0, [s^, • However there is no 
direct evidence to confirm that this particular pathway 
is that taken either by biological or laboratory-prepared 
model membranes. 

In light of the above, it would certainly be desirable 
to examine the fusion pathway in a system whose com- 
ponents are described by a microscopic model. Such 
examination has begun recently. A minimal model, 
consisting of rigid amphiphiles of one hydrophilic and 
two hydrophobic segments, and with no explicit solvent, 
was studied with Brownian dynamics simulations |37j . 
At the same time, a model of more complex, flexible, 
chain molecules, widely employed in the polymer com- 
munity, was used by us to study bilayers composed of 
amp hiphilic, diblock copolymers in a hydrophilic solvent 
|34| . Such copolymers are known, in fact, to form bilayer 
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vesicles which can undergo fusion |9j. This system was 
studied by means of Monte Carlo simulations. Both theo- 
retical studies observed the formation of the initial stalk, 
but found that the subsequent fusion pathway was not 
the usual hcmifusion mechanism, but involved interme- 
diates that broke the axial symmetry. In particular, off to 
the side of the initial stalk, the formation of small pores 
in each of the fusing bilayers was clearly seen. (Wc shall 
refer to these small pores which span one bilayer only 
as "holes" to avoid confusion with fusion pores which 
span both bilayers.) It is intriguing that the two stud- 
ies observed the same fusion pathway even though the 
architecture of the constituents of the two systems dif- 
fered considerably, sharing little other than the generic 
property of being amphiphilic and capable of bilayer self 
assembly. 

The two investigations gave a first glimpse of a fusion 
pathway which differs from the hcmifusion mechanism, 
but did not provide a great deal of quantitative detail. 
In this paper we present an extensive study of the same 
microscopic model we employed previously, and offer suf- 
ficient quantitative evidence to substantiate our earlier 
observations. 

Naturally we are concerned with the question of 
whether the fusion pathway we observe in our model sys- 
tem is relevant to membrane fusion in biological systems. 
The architecture of the components in our system obvi- 
ously differs greatly from those of biological lipids, and 
it is not clear how one should compare the systems. We 
make such an attempt by calculating several dimension- 
less ratios which can be formed from membrane parame- 
ters and comparing those in our system with ratios char- 
acteristic of vesicles formed of block copolymers, and of 
liposomes. (See Table I below.) Ultimately we cannot be 
sure of the systems to which our results apply and under 
what conditions, save the very particular ones that we 
have simulated for the particular case of block copoly- 
mers. In this sense, our results must be evaluated in 
the same way as those from the phenomenological the- 
ories; they must be compared to experiment. We do so 
in the Discussion. In particular we note that our mech- 
anism predicts that the fusion rate depends on lipid ar- 
chitecture and membrane tension, that there is mixing of 
lipids in the cis leaves before mixing of contents, and that 
there is also mixing of lipids between cis and trans layers. 
Of most interest, our mechanism predicts that transient 
leakage is causally linked to the process of membrane fu- 
sion. 



II. SIMULATION DETAILS 

Simulation of membrane fusion in a fully chemically 
realistic model would be most valuable, because it could 
provide information about specific structural changes on 
the atomic level. This would be particularly important 
if changes in molecular conformations entailed a qualita- 
tive spatial redistribution of hydrophilic and hydrophobic 



segments. Unfortunately, the simulation of atomistically 
faithful models can only follow the time evolution of a 
few hundreds of lipid molecules over a few nanoseconds 
even on state-of-the-art supercomputers. Given that the 
time scale of membrane fusion is on the order of millisec- 
onds and involves lengths on the order of a few tens of 
nanometers, an atomistic simulation of the fusion process 
is not yet feasible and one has to resort to coarse-grained 
models. 

Coarse-grained models of amphiphilic chain molecules 
have been used with great success to investigate common 
features of self-assembly. Such models retain only those 
molecular properties that arc necessary for self-assembly, 
such as the connectivity of hydrophilic and hydrophobic 
portions along the amphiphilic molecule, and the mu- 
tual repulsion between these different kinds of segments, 
and ignore specific chemical or electrostatic interactions. 
The usefulness of this approach rests on the observation 
that chemically very different systems, such as biological 
lipids in aqueous environment and block copolymers in 
a homopolymer environment, exhibit a common phase 
behavior and similar structural patterns on length scales 
comparable or larger, than the molecular size. The self- 
assembly of amphiphiles into bilayer membranes itself is 
an example of a universal behavior, i.e. one which does 
not depend on fine details of the underlying architecture. 
It has been successfully studied by coarse-grained mod- 
els We expect that all membranes can be caused 
to fuse, however there may be several different pathways 
which are taken by different systems under different con- 
ditions. Our purpose here is to demonstrate one path 
which is taken in a system modeled microscopically. 

We employ the bond fluctuation model of a polymer 
chain, which has been used previously to study pore nu- 
cleation in a symmetric bilayer membrane under tension 
[33. Much is known about the structure and thermody- 
namics of this model, and the parameters can be mapped 
onto the standard Gaussian chain model of a dense mix- 
ture of extended molecules. In this three-dimensional lat- 
tice model, each segment occupies a lattice cube. No two 
occupied cubes can share any corner, a rule that mim- 
ics hard-core repulsion interaction. Furthermore this en- 
sures that the lattice spacing is sufficiently smaller than 
the width of interfaces so that the effect of the discretiza- 
tion of space is minimal. To ensure that the chain of 
segments can not intersect itself, the segments are con- 
nected by bond lengths that cannot be too large. In par- 
ticular, neighboring segments along the chain can be con- 
nected by one of 108 bond vectors of lengths 2, -s/S, V6, 3 
or -s/TO measured in units of the lattice spacing u. The 
angles between adjacent monomers can take on any of 
87 values. The large number of bond vectors and the 
extended segment shape allow a rather faithful approxi- 
mation of continuous space, while retaining the compu- 
tational advantages of lattice models. The amphiphilic 
molecules consist of 11 hydrophilic segments and 21 hy- 
drophobic segments. This asymmetry mimics the ratio of 
head and tail size in biologically relevant lipid molecules. 
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and is slightly smaller than employed by us previously 
[sJl- We reduced, in this study, the asymmetry of the 
molecules so that a solvent-free system not only would 
be in a lamellar phase (La), but would also be further 
than in our previous study from the boundary separating 
the lamellar and inverted-hexagonal (Hn) phases. The 
solvent in our system is represented by chains of 32 hy- 
drophilic segments, i.e. we conceive a hydrophilic chain 
as a small cluster of solvent molecules, just as in other 
coarse-grained modeling ji^l The mean head-to-tail dis- 
tance of the amphiphiles and solvent molecules is I7u. 
Like segments attract each other and unlike segments 
repel each other via a square well potential which com- 
prises the nearest 54 lattice sites. Each contact changes 
the energy by an amount e = O.lTGSOfcsT. This corre- 
sponds to an intermediate segregation xiV w 30 in terms 
of the Flory-Huggins parameter x- If wc increased the 
incompatibility much more, we would reduce the inter- 
facial width between hydrophilic and hydrophobic seg- 
ments to the order of the lattice spacing and the local 
structure of the lattice model would become important. 
If we decreased the incompatibility, we would reduce the 
clear segregation between hydrophilic and hydrophobic 
regions. Similarly if we replaced the solvent homopoly- 
mers by monomers, we would effectively reduce the in- 
compatibility |28|, and again reduce the segregation be- 
tween the diblock and solvent hydrophilic segments and 
the diblock hydrophobic segments. Were we to increase 
the incompatibility to restore the desired degree of segre- 
gation, we would again reduce the interfacial width of the 
membrane to an extent that lattice effects would become 
important. 

Monte Carlo simulations are performed in the canon- 
ical ensemble, except for some runs described in section 
III. The segment number density, i.e. the fraction of lat- 
tice cubes occupied by segments, is fixed at p = 1/16. 
The conformations are updated by local segments dis- 
placements and slithering-snake-like movements. The 
different moves are applied with a ratio 1:3. We count 
one attempted local displacement per segment and three 
slithcring-snakc-like attempts per molecule as four Monte 
Carlo steps (MCS). This scheme relaxes the molecular 
conformation rather efficiently. The latter moves do not 
mimic a realistic dynamics of lipid molecules and we 
cannot identify straightforwardly the number of Monte 
Carlo steps with time. The density of hydrophilic and 
hydrophobic segments, however, is conserved so that the 
molecules move diffusively. Moreover, the molecules can- 
not cross each other during their diffusive motion. In 
this sense we have a slightly more realistic time evolu- 
tion on local length scales than in dissipativc particle 
dynamics simulations |43l |. but Monte Carlo simulations 
cannot include hydrodynamic flows, which might become 
important on large length scales. At any rate, we do not 
expect the time sequence to differ qualitatively from that 
of a simulation with a more realistic dynamics on time 
scales much larger than a single Monte Carlo step. Most 
importantly, fusion is thought to be an activated process, 



therefore the details of the dynamics only set the abso- 
lute time scale, but the rate of fusion is dominated by 
free energy barriers encountered along the fusion path- 
way, which are independent of the actual dynamics used. 



III. PREPARATION AND PROPERTIES OF A 
SINGLE BILAYER 

It seems clear that bilayers that are under no stress 
will not undergo fusion, as there is no free energy to be 
gained by doing so. So to promote fusion, we have sub- 
jected the studied bilayers to lateral tension. This has 
been done by providing the system with fewer molecules 
than arc needed to span the given area of our sample 
cell with bilayers that are tensionless. Of course we need 
to know just how many molecules are needed to make 
a tensionless bilayer that spans the cell. To determine 
properties of the tensionless bilayer, we made use of the 
definition of the tension in this liquid-like bilayer as the 
derivative of the free energy with respect to the bilayer 
area at constant temperature and particle number. We 
therefore investigated an isolated bilayer with a straight, 
free edge. A simulation cell of size 64u x 200u x 64u with 
periodic boundary conditions in all dimensions was used. 
The bilayer, oriented in the x — y plane, spanned the sys- 
tem in the short, x, direction, but did not span the system 
in the long, y, direction. Its extension in this direction 
adjusted itself until it neither grew nor shrank. Thus the 
surface tension, 7, of the bilayer was zero. This vanish- 
ing value includes, of course, the contributions from the 
fluctuations of the bilayer. Even though the tension van- 
ished, these fluctuations of the mid-plane were not very 
large due to the stiffness of the rather small patch of 
membrane considered. A typical snapshot of the bilayer 
configuration is shown in Fig. ^(s)- Rearrangement of 
amphiphiles at the bilayer free edges is clear. The aver- 
age profile along the y axis, the long axis of the bilayer, 
is presented in panel (b). To obtain it, we have averaged 
the profiles along the x and z direction and estimated the 
instantaneous angle the bilayer makes with the z direc- 
tion (to correct for the difference between projected and 
true area). We observe for these laterally averaged pro- 
files that the edge of the bilayer is slightly thicker than 
the middle, the increase is about 7% for the amphiphilic 
segment density and about 16% for the density of hy- 
drophobic segments. Away from the edge, the densities 
decay exponentially to those of the uniform bilayer {i.e 
without an edge), and we estimate the thickness of the 
tensionless bilayer from that in the middle, finding it to 
be do — 31w. 

The profiles across a single bilayer of thickness do = 31 
are shown in Fig. H^c). They were obtained by simula- 
tion in a cell 40u x 40u x 80u in which the bilayer spanned 
both short directions. One sees that hydrophobic and hy- 
drophilic regions are clearly separated, but there is some 
intcrdigitation of the hydrophobic tails emerging from 
the opposing monolayers. 
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Knowing the thickness of the tensionless bilayer, we 
know the number of molecules needed to span the sim- 
ulation cell with such a bilayer, and can control tension 
by varying the number of molecules introduced into the 
cell. We cannot determine this tension, as one might 
in a Molecular Dynamics simulation, from the excess 
tangential stress in the intcrfacial zone because we em- 
ploy a lattice model. Nevertheless, we can determine 
the tension purely from thermodynamic relations. To 
do so, we assembled a single bilayer in a system of size 
156it X 156u X 64u, where the bilayer spanned the system 
in the x — y plane. Using semi-grand canonical iden- 
tity switches between amphiphiles and solvent, we con- 
trolled the exchange potential A/x between the species 
and monitored the thickness of the bilayer (measured 
by the areal density of amphiphiles). The dependence 
of the bilayer tension on the chemical potentials of the 
amphiphile, /ic, and solvent, is given by the Gibbs 
absorption isotherm |^: 



L^d7 = —Sncdjic — SuAdiiA ~ —SncdA^ (1) 



where Snc and Sua are the excess number of molecules in 
the bilayer. In the last approximation we have assumed 
that the liquid is incompressible i.e. Sua ~ —Snc, and 
the solubility of the amphiphile in the solvent is vanish- 
ingly small. Results of the simulation for the number of 
amphiphiles dnc as a function of the exchange poten- 
tial A/i ~ He — fJ-A are shown in Fig. ^d) . Using the 
thickness of the tensionless bilayer, we can estimate the 
tension of an arbitrary bilayer as a function of exchange 
potential or of thickness by integrating Eq. (Q. The 
results, in reduced units of bare A-B homopolymer inter- 
facial tension 70 = O.OGSfcsT/w^, are shown in the inset 
of Fig. Ekd). Dashed lines in Fig. ^d) and the inset 
correspond to the tensionless bilayer. Comparison of the 
relevant structural and elastic properties of the polymer- 
somes, liposomes and simulated membranes is provided 
in Tabled 

We are now in a position to simulate bilayers under 
a given tension in the canonical ensemble. Knowing the 
area of our simulation cell, and the segment density, we 
add the number of amphiphiles which will produce a bi- 
layer of a given thickness. From Fig. Q^d), wc know 
what tension is placed on this bilayer. For our study of 
two bilayers under tension, we have chosen their thick- 
ness to be c? = 25u, smaller than the thickness do = 3lu 
of the tensionless bilayer. This corresponds to a ten- 
sion of the order of 7/70 ~ 0.75 and an area expansion, 
AA/Aq « 0.19. We know from our simulations that a sin- 
gle bilayer of the thickness chosen, d = 25u, is metastable 
on the time scale of fusion, i.e. the small holes, which 
appear transiently, do not grow past their critical size on 
the time scale of fusion in our simulations. 



IV. PREPARATION OF A SYSTEM OF TWO 
BILAYERS 



We begin with a system containing only amphiphiles. 
It is 156u X 156w x 25u with periodic boundary conditions 
in the long directions, and hard, impenetrable, walls in 
the short direction. They attract the hydrophilic por- 
tion of the amphiphile and repel the hydrophobic portion. 
These interactions extend over two layers nearest to the 
wall, and each contact changes the energy by Q.GksT. 
The amphiphiles assemble into a bilayer structure which 
is free of defects. 

Two of these flat bilayers are then stacked on top of 
each other with a distance of A between them, and are 
embedded into a simulation cell with geometry 156w x 
156m X 126m. There are no walls at this point, and pe- 
riodic boundary conditions are utilized in all three di- 
mensions. The conditions of flat bilayers mimic the ap- 
proach of two vesicles whose radii of curvature are much 
larger than the patch of membrane needed for fusion. 
The solvent of homopolymers is then inserted into the 
simulation cell via grand canonical, configurational bias, 
Monte Carlo moves at infinitely large chemical potential 
of the homopolymer until the segment number density of 
p = 1/16 is reached. The initial distance A between the 
bilayers translates into the thickness of the residual sol- 
vent layer between the two membranes. We have carried 
out the most extensive series of runs with A = 10m and 
unless specified otherwise, all our results are for that sep- 
aration. Because the solvent homopolymers are flexible 
coils, and each represents a cluster of solvent molecules, 
many layers of solvent segments are represented between 
the bilayers at this separation. In our previous simu- 
lations |3J|, we set A = and observed qualitatively 
similar behavior as we do with this larger separation. 
We increased the separation for this extensive study be- 
cause, as expected, the rate of fusion events decreased 
(see Sec. V. below) and this allowed us to observe the 
sequence of structural rearrangements more clearly than 
in our previous work. The separation chosen, a bit less 
than half the thickness of one bilayer, is comparable to 
the separation at which fusion occurs when mediated by 
hemagglutinin [13. A snapshot of the two bilayers is 
shown in Fig. [3 Hydrophobic and hydrophilic segments 
of amphiphiles are shown as dark and light gray spheres. 
For clarity solvent segments, which are present in the 
simulation, are not shown 

Thirty-two independent starting configurations were 
prepared, each containing 194, 688 segments correspond- 
ing to about 3,613 amphiphiles and 3,708 solvent 
molecules. After every 25,000 Monte Carlo steps, a con- 
figuration was stored for further analysis. Ten thousand 
hours of CPU time were utilized in the course of this in- 
vestigation, with thirty two processors running for about 
two weeks. 
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THE PROCESS OF FUSION 



The stalk and associated hole formation 



It is straightforward to monitor the internal energy of 
the system during the simulation because this energy 
arises solely from contacts between segments, and the 
locations of all segments are known. (In contrast, the 
free energy cannot be obtained directly.) We show in 
Fig.|21the behavior of the internal energy of two systems, 
one separated by a distance A = 4m (squares), and the 
other with A = lOu (circles). The energy is plotted in 
units of fc^T, as a function of time, in units of 25,000 
Monte Carlo steps. The energy initially decays, which 
reflects the equilibration of the system. During this ini- 
tial relaxation of the starting configuration the interface 
between the bilayer and the solvent adjusts locally. The 
time scale of this initial relaxation (less than 25, OOOMCS) 
is independent of the distance between the bilayers, and 
is about two orders of magnitude smaller than the time 
scale on which the fusion pore forms. Due to this separa- 
tion of time scales between initial relaxation and fusion 
we do not expect the preparation of the starting con- 
figuration to affect the fusion process. Similarly we do 
not expect our results to depend on our particular choice 
of relaxation moves, as other choices would also lead to 
relaxation of the bilayers which takes place on a much 
shorter time scale than does fusion. 



After the initial relaxation, two subsequent time 
regimes can be identified. First the energy rises slowly. 
Two mechanisms contribute to this increase of the en- 
ergy. On the one hand, capillary waves of the hy- 
drophilic/hydrophobic interfaces become thermally ex- 
cited. They increase the effective interface area and 
thereby lead to a slow increase of the energy. Addition- 
ally, undulations result in the formation of stalks and 
holes. We shall discuss the details of this process in the 
next subsection. Later, around 320 x 25000 MCS the en- 
ergy decreases rapidly. This final decrease of the energy 
results from the fusion of the membranes which releases 
some of the tension stored in them. As noted above, the 
fusion occurs more rapidly the closer the bilayers, as ex- 
pected. The increase in energy preceding fusion reflects 
the formation of fusion intermediates, the focus of our 
study. 

The inset of Fig. |21 shows the fluctuations in the en- 
ergy, i.e., the fluctuations between the 32 different runs 
at equal time. Strong fluctuations indicate energy differ- 
ences between the independent runs. The peak at around 
t ^ 320 indicates that some systems have already formed 
a fusion pore (and therefore have a lower energy) while 
others systems have only stalks and holes (and therefore 
have a higher energy). The vertical bar indicates the time 
we have chosen to indicate on several figures the onset of 
fusion. The width of the peak provides an estimate for 
the spread of the time at which a fusion pore appears. 



During the initial stage of simulations the fluctuating 
bilayers collide with one another frequently and some- 
times form small local interconnections. For the most 
part, these contacts are fleeting. Occasionally we ob- 
serve sufficient rearrangement of the amphiphiles in each 
bilayer to form a configuration, the stalk, which connects 
the two bilayers, (see Fig.^Ja)) and which is not as tran- 
sient. Such a stalk was hypothesized l ong ago to be in- 
volved in the initial stages of fusion |l7l |27| . In con- 
trast to stable arrays of stalks which have been observed 
in block copolymer melts ^3 and in lipid systems [48l |. 
those we see are isolated, and increase the free energy of 
the system. We infer the latter from two observations: 
that the appearance of stalks is correlated with the in- 
crease in the internal energy of the system as a function 
of time shown in Fig. |31 that some stalks vanish with- 
out proceeding further to a fusion pore. Thus it appears 
that the stalk represents a local minimum along the fu- 
sion pathway. Density profiles of the hydrophilic and hy- 
drophobic parts of the amphiphiles in the presence of the 
stalk, and obtained by averaging over configurations, are 
shown in Fig. [S| The dimples in the membranes at each 
end of the stalk axis are notable. What can barely be 
seen is a slight thinning of each bilayer a short distance 
from the axis of the stalk. 

After stalks are formed, the rate of formation of holes 
in either of the two bilayers goes up markedly. This can 
be seen in Fig. |B| in which we plot the fractional area of 
holes as a function of time for the system of two apposed 
bilayers, and compare it to the rate of hole formation in 
an isolated bilayer. In contrast to the large increase in 
the area of holes formed in the apposed bilayers at "time" 
t = 200 X 25000 MCS, the fractional area in single bilay- 
ers fluctuates somewhat about an average value which 
is rather constant over time at a value of approximately 
0.004. Comparison with Fig. 13 shows that the increase 
in the rate of hole production in the apposed bilayers in 
this system with bilayer spacing A = lOu is correlated 
in time with the decrease in the energy of the system, 
and it is reasonable to infer that the decrease in energy 
is caused by the production of holes and, later, the fusion 
pore. Similarly, during the time before this increase in 
hole production, stalks are forming, and it is also reason- 
able to infer that the increase in energy is due to their 
formation. 

The locations of stalks and holes are correlated; holes 
form close to the stalks, and the stalk elongates and 
moves to surround the hole. A snapshot of this is shown 
in Fig. |7Ja) and (d). In both snapshots an elongated 
stalk is seen and a small hole is formed in the upper bi- 
layer next to the stalk. The extent to which holes are, on 
average, found close to a stalk can be determined from 
the holc-stalk correlation function 

T,r,.rJi\^s-rh\-r)Pshirs,rh) 
E,....^(|r.-r,.|-.) 
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where Pshi^^s, r/i) is the joint probabihty that the lateral 
position is part of a stalk and is part of a hole, and 
S{r) is the Dirac delta function. The value of g(r) at large 
distances r is proportional to the product of the areal 
fraction of holes and stalks. This correlation function 
is shown in Fig.|Hl The scale of g(r) increases with time 
indicating the simultaneous formation of stalks and holes. 
The figure shows that the correlation peaks at a distance 
of about 16u, and falls rapidly at larger distances. (Recall 
that each bilayer has an average thickness of 25m.) 

It is not difficult to understand why the presence of a 
stalk promotes hole formation. First, if the hole forms 
close to a stalk, then the line tension, or energy per unit 
length A, of that part of the hole near the stalk is signif- 
icantly reduced. This can be seen from the schematic in 
Fig. |5| In the upper part of the figure, we show a hole 
which has formed far from a stalk, while in the lower, 
we show a hole which has formed close to one. It seems 
clear that the line tension in the latter is reduced simply 
due to the reduction of curvature of the hydrophobic- 
hydrophilic interface. The second reason that the stalk 
formation encourages the appearance of holes is due to 
the slight thinning of the membrane in the vicinity of the 
stalk to which we alluded earlier. Further it has been 
suggested recently that the local surface tension in the 
neighborhood of a defect, such as a stalk, is increased 
significantly making such a location the likely site of 
hole formation. 

Now that one hole has formed next to the stalk, and 
the stalk has begun to surround it, two other events occur 
to complete the formation of the fusion pore. They are: 
i) a second hole forms in the other bilayer, ii) the stalk 
surrounds the hole(s) to form the rim of the fusion pore. 
We have observed these steps to occur in either order, 
and will briefly discuss them separately. 



B. Pathway 1. Rim formation followed by 
appearance of a second hole 

In this scenario, a hole appears in one bilayer and the 
stalk completely surrounds it rather rapidly. A snapshot 
of the system in this configuration is shown in Fig. [7Kb) . 
This looks very much like a hemifusion diaphragm which 
has been suggested by many authors as an intermediate 
stage in fusion 0, 123, 0|. However, this diaphragm is 
quite different from the usual hemifusion one that con- 
sists of two trans monolayers of the fusing membranes. 
In contrast, the diaphragm we observe is made of one 
of the pre-existing bilayers; that is, it is made of cis and 
trans leaves. The appearance of a hole in this diaphragm, 
as shown in Fig. EKc), and its expansion completes the 
formation of the fusion pore. 



C. Pathway 2. Appearance of second hole followed 
by rim formation 

In this scenario, a hole appears in one bilayer and, 
before the stalk completely surrounds it, a second hole 
appears in the other bilayer. The stalk tries to surround 
them both, and aligns them in doing so. In Fig. [3(e) 
we show one stage in this process. One sees a large hole 
in the upper bilayer. A small hole is formed in the lower 
bilayer next to the stalk. Eventually, the stalk aligns and 
completely encircles the holes (sec Fig.[7Kf)) to form the 
final fusion pore shown in Fig. ^b). Again, the driv- 
ing force for the stalk to surround the two holes is the 
reduction in their (bare) line tension. Because the stalk 
aligns and surrounds two holes, we observe this pathway 
to be somewhat slower than that of pathway 1 in which 
the stalk need only surround one hole. 

Once the fusion pore has formed, by either of the above 
mechanisms, it expands, driven by the reduction in sur- 
face tension. The growth of the fusion pore eventually 
slows and ends as the pore reaches its optimum size de- 
termined by the finite size of our cell. 



VI. DISCUSSION 

We have carried out extensive Monte Carlo simula- 
tions on the fusion of two bilayer membranes comprised 
of amphiphilic molecules immersed in solvent. The am- 
phiphiles and solvent are modeled by copolymers and ho- 
mopolymers respectively. The membranes are under ten- 
sion. The mechanism of fusion that we see begins with 
a stalk, as posited years ago, and incorporated in almost 
all fusion scenarios. However what follows after stalk 
formation is different from all other mechanisms which 
have been proposed save that presented independently 
by Noguchi and Takasu [s^- In particular, the fusion in- 
termediates we see break the axial symmetry which has 
been assumed in almost all previous calculations. We 
observed that the stalk destabilizes the bilayers by cat- 
alyzing the creation of small holes in them. We argued 
that the mechanism behind this is quite simple: the en- 
ergy per unit length of the edge of a hole is reduced when 
the edge is adjacent to a stalk. For the same reason, the 
stalk will try to surround the hole formed in one bilayer 
once the hole has appeared. Two slightly different path- 
ways to the final fusion pore were observed differing only 
on whether the hole in the second bilayer, which is neces- 
sary for complete fusion, appears before or after the stalk 
completely surrounds the first hole. 

The question now arises as to whether the pathway we 
see in the model system is that which occurs in biologi- 
cal fusion. There are many differences between the model 
studied and a biological system. Perhaps the most obvi- 
ous is that we have modeled fiexible, single chain block 
copolymers, not lipids with two semiflexible tails and a 
rigid head. How is one to determine whether these archi- 
tectural differences are significant? It is useful to recall 
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that phcnomenological theories completely ignore the ar- 
chitecture of the membrane constituents and encapsulate 
their effects in a small number of parameters which enter 
the theory, such as the monolayer spontaneous curvature 
and bending modulus. In that same spirit, we can extract 
from our simulation those same parameters and compare 
dimensionless ratios of them to those of other systems. 
We have done that, and presented the results in Table I. 
One sees that the values we obtain are reasonable. The 
ratio of the bilayer compressibility modulus to the hy- 
drophilic/hydrophobic interface tension, Ka/^o^ closer to 
that of liposomes than of polymersomes. The reverse 
is true for the ratio of the monolayer bending modulus 
to the product of surface tension and the square of the 
hydrophobic thickness Kh/jod'^. One line in the table de- 
serves comment, that for the experimental values of the 
bilayer area expansion, AA/Aq quoted at rupture, (the 
"critical values"). That for liposomes is smaller than that 
for polymersomes at rupture, which is equal to the bilayer 
area expansion we utilized. However the values quoted 
at rupture have no thermodynamic meaning, because any 
membrane under tension is inherently unstable and will 
be observed to rupture if the time scale of observation is 
sufficiently long. The experimental values quoted apply 
over some, unspecified, laboratory time scale. On this 
point we add that, as in experiment, we found many of 
our bilayers to rupture over the time we observed them, 
but the time scale for this to happen was significantly 
greater than that for fusion. If the bilayer area expan- 
sion, or equivalently, its tension were reduced, either in 
experiment or in our simulation, the time scales for the 
bilayers to fuse and later to rupture would both increase, 
perhaps to the extent of making impossible the observa- 
tion of fusion. Indeed we chose the value of tension in 
the simulation such that fusion could be observed con- 
veniently. One could still ask whether, in addition to 
increasing the time scale for fusion, a significant reduc- 
tion in bilayer tension would favor an alternative fusion 
pathway. To attempt to answer this question, one could 
contemplate even longer Monte Carlo runs on membranes 
under less tension. 

There are other physical parameters which might af- 
fect the fusion pathway but which are not encompassed 
by the quantities in Table I. For example, one might ask 
whether the fusion pathway is expected to be the same 
for large virus-encapsulating endosomes as it is for small 
synaptic vesicles. Thus one would consider the dimen- 
sionless ratio of the membrane's hydrophobic thickness 
to the radius of the vesicle in question. We have con- 
sidered the simplest case of planar membranes for which 
this ratio is zero. For endosomes encapsulating influenza 
viruses with an average diameter of lOOnm, the ratio is 
small, less than, but of the order of, 0.03, but for synaptic 
vesicles of typical diameter 50nm, it is at least twice this. 
It is not difficult to imagine that for a sufficiently large 
value of this ratio, which implies a small area of contact 
between the fusing vesicles, there might be insufficient 
room for the growth and movement of the stalk we have 



observed, so that our mechanism would be supplanted by 
another. But we do not know this. 

Ultimately the most meaningful test of the applicabil- 
ity of our mechanism to biological fusion is comparison 
to experiment, and our scenario does have testable con- 
sequences. First, because of the initial stalk formation, 
one expects to see the mixing of lipids in the two proxi- 
mal layers before the fusion pore opens, if it forms at all, 
a result which is in accord with experiment ^3 1^ ■ 
Second, due to the formation of holes in each bilayer 
near a stalk, our scenario allows for the mixing of those 
lipids in the cis and trans leaves of one bilayer and also of 
lipids in the cis leaf of one bilayer with those in the trans 
leaf of the other. The standard hemifusion mechanism 
does not permit either process. Note that this move- 
ment is different from lipid flip-flop which is known to 
be very slow. Mixing of lipids between the cis and trans 
monolayers has been observed in fusion p^ . , but it 
has not yet been determined from which membrane they 
originate and in which membrane they terminate. We 
have monitored the amphiphiles to see whether they re- 
main in the leaf in which they were situated at the be- 
ginning of the Monte Carlo run, or mix with amphiphiles 
in other leaves. Instantaneous assignment of amphiphiles 
to a respective monolayer was determined by the center 
of mass of their hydrophilic part. The results are shown 
in Fig. ^1 They share with experiment the fact that 
the membrane of origin is not distinguished nor is the 
membrane of final residence. In order to evaluate the 
results for the apposed bilayers under tension, we have 
also included those for the single isolated bilayer under 
zero tension and under the same tension (7/70 = 0.75) 
as in the simulations of fusion. Lateral tension greatly 
enhances the flip-flop rate in the single bilayer system. 
This effect can be explained by an overall thinning of 
the membrane, which lowers the translocation barrier, as 
well as by the diffusion of amphiphiles through the tran- 
sient holes formed under tension. In the simulations of 
the apposed bilayers, translocation of amphiphiles from 
the trans leaves initially follows the same dynamics as in 
the single bilayer system, but eventually deviates from 
it, apparently due to the formation of holes facilitated 
by the appearance of the stalks, as discussed in the pre- 
vious section. Amphiphiles from the cis leaves undergo 
mixing to the largest extent, as would be expected due to 
stalk formation. Third, our mechanism allows for tran- 
sient leakage during fusion. As noted earlier, there will 
be greater leakage if fusion occurs via pathway two in 
which the stalk aligns and surrounds two holes than if it 
occurs via pathway one in which the stalk rapidly sur- 
rounds one hole before the second appears. Clearly the 
amount of leakage depends on the size of the transient 
holes formed in the bilayer, the time between the forma- 
tion of the initial stalk and the completion of the fusion 
pore, and the diffusion constant of the molecules which 
leak. This constant introduces another time scale whose 
magnitude, relative to that of fusion pore formation, de- 
termines whether the fusion process is observed to be 
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leaky or tight. 

It is clear that within our mechanism, leakage via tran- 
sient holes and fusion via pore formation are correlated 
in space and time. The latter is shown in Fig. II II which 
presents, as a function of time, the area of holes and that 
of fusion pores from one of the simulation runs. One 
sees in this figure, as in the Monte Carlo snapshots, that 
the rate at which holes appear, and therefore the rate at 
which leakage should occur, increases significantly before, 
and is correlated with, the formation of fusion pores. . 
Once the fusion pore has formed, the creation of other 
holes decreases due to release of tension initially stored 
in the membranes. 

The question of whether transient leakage is charac- 
teristic of membrane fusion is an open one. On the one 
hand, some experiments detect no leakage |^ l47l |. 
while on the other there is a great deal of evidence that 
fusion of biological membranes is, indeed, a leaky pro- 
cess d EH EE EH Ell- It could be argued that observed 
leakage is due to the presence, in these experiments, of fu- 
sion proteins, such as hemagglutinin, which arc certainly 
present in the vicinity of fusion, and which are known 
to undergo conformational changes in which part of the 
protein inserts itself into the target vesicle. In support of 
this view, one could cite the well known ability of fusion 
peptides to initiate erythrocyte hemolysis [s^. Such pep- 
tides are not included in our model. This argument is 
vitiated, however, by the observation that leakage is also 
detected in the fusion of model membranes without such 
peptides 0,0,131 ■ these experiments, large molecules 
such as polyethyleneglycol, are used to bring the fusing 
vesicles together. It would be difficult to argue that these 
molecules, which undergo no conformational change, are 
responsible for the leakage as they generate an attractive 
osmotic force between the vesicles precisely because their 
large size makes it difficult for them to enter the region 
where the vesicles are closely opposed. 

One test that might distinguish whether leakage simply 
accompanies fusion or is causally related to it is provided 
by the observation above that in our mechanism tran- 
sient leakage is correlated in time and space with fusion. 
Just such an experiment to measure these correlations 
has been carried out recently E3 , is reported in the 
paper accompanying this manuscript. They observe that 
leakage is, in fact, correlated spatially and temporally 
with the process of fusion. Indeed, their results com- 
paring the time sequence of the electrical conductance 
arising from leakage with that arising from fusion, shown 
in their Fig. 5 displays a remarkable similarity to our 
results comparing the time sequence of the areal fraction 
taken up by holes with that taken up by fusion pores, 
our Fig. 11. 

While the congruence between the predictions of our 
model and experiment are very encouraging, there are 
further tests we should like to apply to it. Foremost 
among these is to determine the free energy barriers for 
the various steps along the fusion pathway. As noted 
above, it is relatively simple to determine the internal en- 



ergy during the course of the simulation as one need only 
monitor the interactions between all segments. But the 
simulations cannot easily evaluate the entropy changes 
along the fusion pathway or, therefore, the free energy 
barrier. To determine the actual value for the free en- 
ergy barrier, calculations using self-consistent field the- 
ory, which have been extremely successful in describing 
the phase behavior of amphiphiles j29|, i3(3 | are currently 
being pursued by us. Also, elastic constants of the simu- 
lated amphiphilic monolayers, e.g. calculated inlsE could 
be employed in the simpler phenomenological theories, 
which have proved to be so useful. Comparison with the 
full self-consistent field calculations would permit deter- 
mination of the accuracy of these elastic models in de- 
scribing the highly curved intermediates involved in the 
fusion reaction. Furthermore, there is an extensive ex- 
perimental evidence on the effect of lipids of differing 
architecture on fusion [E E3- The self-consistent field 
theory is able to describe such differences 25, 28J and to 
determine both the spatial distribution of different am- 
phiphiles in inhomogeneous structures such as the stalk, 
the holes and the fusion pore, as well as the change in 
the free energy of these structures. Results of these in- 
vestigations will be published separately. 

It would be of great interest to repeat our simulations 
under different membrane tension, as this would help 
to clarify the importance of fusion peptides in bring- 
ing about such tension. Finally, it would be desirable 
to carry out simulations in which fusion peptides are 
included explicitly. One could investigate whether the 
membrane perturbations associated with such inclusions 
provide sites for the nucleation of the small holes that 
are necessary for the formation of the fusion pore. If 
this were so, one could test the further inference that, by 
providing nucleation sites in close proximity, one in each 
membrane, such peptides facilitate sucessful and rapid 
fusion thereby reducing leakage. 
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TABLE I: Structural and elastic properties of bilayer mem- 
branes: dc - thickness of membrane hydrophobic core, / - 
hydrophilic fraction, Co - monolayer spontaneous curvature, 
AA/Ao - bilayer area expansion (critical value for the ex- 
perimental systems, and the actual strain used in simula- 
tions), Ka - bilayer area compressibility modulus, Kb - mono- 
layer bending modulus, 70 - hydrophilic/hydrophobic inter- 
face tension (oil/water tension of 50pN/nm for the experi- 
mental systems, and A/B homopolymer tension for the sim- 
ulations). Data on E07 po lym ersomes is taken fromj^ and 
on lipids from (a): El (b): E3, (c): H and (d): ^ (see also 
^ttp://aqueous.labs.brocku.ca/lipid/ ). Values of dc, Co and 
Ka for DOPE were obtained by linear extrapolation from the 
results on DOPE/DOPC(3:l) mixtures and pure DOPC. Val- 
ues of Kb, 70, and Co for the simulated model were calculated 
by us using the method of.33t 




FIG. 1: (a) Snapshot of an isolated bilayer in the tensionless state. Hydropliobic and liydrophihc segments of amphiphiles are 
shown as dark and light gray spheres. Solvent segments are not shown for clarity, (b) Density profiles along the y axis. The 
edge of the bilayer is thicker than its middle, (c) Profiles across the bilayer for a lateral patch of size 40m x 40u. (d) Dependence 
of the bilayer thickness on the exchange chemical potential between amphiphiles and solvent. The inset displays the tension 
7 of the bilayer as a function of exchange potential. 




FIG. 2: Snapshot of the initial configuration in the two bilayer system. Hydrophobic and hydrophilic segments of amphiphiles 
are shown as dark and light gray spheres. Solvent segments are not shown for clarity. 
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FIG. 3: Evolution of internal energy in fusion simulations. The two curves correspond to initial bilayer separations A — 
(squares) and A = lOu (circles). To reduce fluctuations, the data are averaged over all 32 configurations at equal time and 
additionally over small time windows. The large negative value of the energy mirrors the attractive interactions in the solvent. 
The inset shows the sample-to-sample energy fluctuations as a function of time. Large fluctuations identify the onset of fusion. 
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FIG. 4: Representative snapshots of (a) the stalk intermediate and (b) the complete fusion pore from one of the simulation 
runs. Each configuration is shown from four different viewpoints. Hydrophobic core is shown as dark gray, the hydrophilic- 
hydrophobic interface (defined as a surface on which densities of hydrophilic and hydrophobic segments are equal) is light gray. 
Hydrophilic segments are not shown for clarity. Top- and bottom- left sub-panels have been generated by cutting the system 
along the middle x — y-plane, the top and bottom halves are viewed in the positive (up) and negative (down) z-direction 
correspondingly. Top- and bottom- right sub-panels are side views with cuts made hy x — z and y — z planes correspondingly. 
Grid spacing is 2Qu. 3D-orientation axis are the same for all snapshots and shown in panel (a). 



15 




-80 -60 -40 -20 20 40 60 80 

x/u 

FIG. 5: Density distribution of segments in the stalk, averaged over all simulation runs. At each point only the majority 
component is shown; solvent as white, hydrophobic and hydrophilic segments of amphiphiles as black and gray respectively. 
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FIG. 6: Area of holes vs. time in the system of two apposed bilayers (gray for one bilayer and black for the other on the bottom 
panel) and in an isolated bilayer (top panel). 
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FIG. 7: Two observed pathways of fusion process. The snapshots are taken from two representative simulation runs. Each 
configuration is numbered by the time (in multiples of 25,000MCS) at which it was observed. See Fig. |1] for explanation of 
the graphics shown. For discussion of the mechanism see text. 
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FIG. 8: The hole-stalk correlation function at early times. 




FIG. 9: Schematic explanation of the line tension reduction near the stalk. 




FIG. 10: Probability of finding an amphiphilic molecule in its original monolayer after time t. The solid and dashed lines refer 
to simulations of a single bilayer under tension 7/70 = and 7/70 = 0.75 respectively. The lines with symbols present the 
results obtained in the simulations of fusing bilayers. Squares and and circles refer to cis and trans monolayers respectively. The 
time period corresponds to the formation of stalks and holes. Error bars show standard deviations obtained from thirty-two 
runs. 
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FIG. 11: Area of pore (symbols) and of holes (lines) vs. time for one simulation run (identical to Fig. EJ. Note the different 
scale for pore and hole areas. 



